Further details on the phase diagram of hard ellipsoids of revolution 
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In recent work we revisited the phase diagram of hard ellipsoids of revolution (spheroids) by means 
of replica exchange Monte Carlo simulations. This was done by setting random initial configurations, 
and allows to confirm the formation of sm2 crystal structures at high densities [Phys. Rev. E 75, 
020402 (2007)] for large anisotropics and stretched-fcc for small anisotropics. In this work we 
employed the same technique but setting the starting cells as sm2 crystal structures having the 
maximum known packing density [Phys. Rev. Lett. 92, 255506 (2004)]. This procedure yields 
a very rich behavior for quasi-spherical oblates and prolates. These systems, from low to high 
pressures, show the following phases: isotropic fluid, plastic solid, stretched-fcc solid, and sm2 solid. 
The first three transitions are first order, whereas the last one is a subtle, probably high order 
transition. This picture is consistent with the fact of having the sm2 structure capable of producing 
the maximally achievable density. 



PACS numbers: 64.30.-t, 64.70.mf, 61.30.Cz 

I. INTRODUCTION 

During quite a long time it was assumed that the max- 
imum achievable density of hard ellipsoids was given by 
the stretched face cubic centered (sfcc), with a volume 
fraction of cp m = n/y/lS ~ 0.7405. However, after the 
observation of certain crystal structures capable of sur- 
passing this threshold [1 j, the high density region of the 
original phase diagram [2j [3] has been significantly mod- 
ified [4-6 . Naturally, the updated phase diagram must 
include the crystal structures which, under large com- 
pressions, produce the maximally achievable densities. In 
particular, Pfleiderer and Schilling [4 found that a fam- 
ily of crystals named sm2 (simple monoclinic with two 
orientations) showed smaller free energies than that of 
the sfcc for ellipsoids having large asymmetries, whereas 
the opposite was found for ellipsoids with small asym- 
metries [5]. On the other hand, the structures given by 
Donev et. al. showing the maximally known achievable 
densities are particular cases of the sm2 family [4j [5] . 

At present, the phase diagram of hard ellipsoids with 
large asymmetry shows consistency between the free en- 
ergy predictions and the structure showing the maximally 
achievable density. From low to high densities, there is an 
isotropic fluid, a nematic fluid, and finally a sm2 crystal 
structure which would yield the densest structure under 
infinite pressure. Conversely, the low asymmetry region 
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still looks incomplete. That is, there is an isotropic fluid, 
a plastic crystal [3] , and a sfcc structure of parallel ellip- 
soids which cannot be compressed to reach the maximum 
packing fraction (parallel ellipsoids cannot exceed the fee 
density limit [7] ) . Consequently, the sfcc structure of par- 
allel ellipsoids must somehow distort to increase its den- 
sity under extreme compression, or suffer another phase 
transition to produce the sm2 structure. 

In recent work we revisited the phase diagram of hard 
ellipsoids of revolution by means of replica exchange 
Monte Carlo simulations [6 j. This was done by setting 
random initial conditions, and the results confirmed the 
spontaneous formation of sm2 crystal structures at high 
densities and relatively large anisotropics, and parallel 
sfcc structures for small anisotropics. In the present work 
we employ the same technique but setting the starting 
cells as perfect sm2 crystal structures having the maxi- 
mum known packing density pQ. This is done with the 
hope that these structures are indeed the ones that reach 
the maximum packing density. If so, the obtained re- 
sults should correspond to equilibrium. Anyway, a direct 
comparison with the pressure-density curves from ran- 
dom initial conditions is possible. If curves match, equi- 
librium would be the case. On the contrary, the method 
will be certainly failing in the sense that the ergodic hy- 
pothesis is violated for the given conditions. 

The setting of perfect sm2 structures as initial condi- 
tions allows us to capture a very high pressure sfcc-sm2 
transition for small anisotropics. In this region, we find 
that the sm2 structure holds only at very high pressures, 
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turning first into a parallel sfcc structure and then into a 
plastic-solid when decreasing pressure and density. Ac- 
cording to our results, obtained for small systems, the 
sfcc-sm2 transition is not first order. The existence of this 
sfcc-sm2 transition occurs only above the plastic solid 
phase, for both, oblates and prolates, whereas the sfcc 
structure region vanishes at large anisotropics. 

The paper is organized in four sections. Following this 
brief introduction, the second section describes the em- 
ployed models and methods. In a third section, the high 
pressure phase transition for the small anisotropic el- 
lipsoids is shown, together with a comparison between 
the results obtained by starting from sm2 structures and 
from random initial configurations. Some remarks and 
conclusions are given in a final section. 

II. MODELS AND METHODS 




FIG. 1. a) The figure illustrates the two- folded structure 
of the sm2 structure, b) Calculation of the distance between 
two equally oriented layers, h, for prolates (oblates result is 
identical) . 



A. Hard ellipsoids model 

In order to perform a simulation of hard ellipsoids, we 
need an efficient way to avoid particle overlapping. For 
this purpose, we use an analytical approach for the exact 
hard ellipsoids contact distance. The expression is based 
on the Berne and Pechukas [8 closest approach distance 
(also called hard Gaussian overlap [9 ), and includes a 
corrective term to improve precision. This term was in- 
troduced by Rickayzen to fix the known T-shape Berne 
and Pechukas mismatch. The Rickayzen-Berne-Pechukas 
(RBP) expression reads [TO] 
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Here, cry and a± are the parallel and perpendicular di- 
ameters with respect to the ellipsoid axis of revolution, 
respectively. We define the aspect ratio of the ellipsoids 
as a = ct||/ct_l, such that a > 1 corresponds to prolates, 
and a < 1 to oblates. Analogously, the maximum as- 
pect ratio is defined as 6 = cF max / cF m i ni where cF max and 
drain are the maximum and minimum axes, respectively. 
\ii and \ij are unit vectors along the axis of revolution 
of each particle, and r is the unit vector along the line 
joining the geometric particle centers. Finally, 7 is in- 
troduced pj] to further approach to the exact Perram 
and Wertheim numerical solution [I2j [13] . 7 values are 
given in reference [II . The average difference between 
the analytical approach and the exact numerical solution 



is always below 0.8% for 0.2 < a < 5, for a collection of 
10 8 random configurations (varying r, u^, and The 
equations of state corresponding to the RBP and PW 
were also compared showing no practical differences for 
oblates with S = 5 [TT1. 



B. Donev-sm2 structure 

Donev et. al. pQ described an arrangement of ellipsoids 
able to surpass the maximal density of the sfcc structure. 
This structure is generated starting from a horizontal 
layer (A) of equally oriented ellipsoids, such that 5 ellip- 
soids are fitted in contact inside a face-centered square 
lattice of side L (see Fig. [T^l) ) . The two axes parallel to 
the horizontal plane are given by cry and <tj_, respectively. 
As we are assuming ellipsoids of revolution, the length of 
the axis perpendicular to this plane will be also given by 
cr±. Within these conditions, the length of the square 
lattice is 
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Then, a second horizontal layer (B) is built on top of 
the first one following the same procedure, but rotated 
by 7r/2 around the perpendicular direction. The sec- 
ond layer is also shifted in such a way that the cen- 
ters of the ellipsoids are located in the holes formed by 
the first layer, each ellipsoid of the second layer touch- 
ing four ellipsoids of the first one. This procedure is re- 
peated successively, leading to a stratification in the form 
ABABAB..., where each layer perfectly fits in the holes 
of the other. In both cases (prolates and oblates), the 
obtained structure is denser than the sfcc structure for 
6 < \/3, and leads to a density maximum when the el- 
lipsoids in the face-centered layers touch six rather than 
four in-plane neighbors, which occurs for 6 = VS pQ. 
Thus, for 6 < y/3, we are taking this structure as the one 
with the highest possible density. 
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FIG. 2. a) Perspective views of unit cells for oblates with a = 
l/y/3 (left), and the one obtained after stretching the particles 
and the unit cell (right), using a scale factor of s = 2.484, to 
produce oblates with cry = 1.198 and a± = 3.594 (a = 1/3). 
b) Initial configurations for 120 oblates with a — l/y/3 (left) 
and a = 1/3 (right). Different colors are used to highlight the 
two different orientations. 



The distance between two equally oriented layers, ft, 
can be analytically determined by considering that the 
ellipsoids of two consecutive layers are in contact. This 
calculation can be simplified to a two-dimensional prob- 
lem in the vertical plane (yz-plane of Fig. [T]): We only 
need to determine the contact point between an ellipse 
of axes (cr_L,(j||) and a circumference of diameter a±. 
Fig. [TJd) illustrates the position and orientation of both 
geometrical objects, which correspond to the shaded el- 
lipsoids of Fig. [IJi). Setting the origin of the Carte- 
sian coordinates at the center of the ellipse, the set 
of points that belong to the ellipse can be written as 
r = (y,z) = (1/2) (cry cos(/>, a± sin0). With this pa- 
rameterization, the unit vector perpendicular to the el- 
lipse at the contact point may be expressed as n = 
(ct_l cos </>, ay sin <j>) / (a]_ cos 2 (j) + a 2 sin 2 (j)) 1 / 2 . The vector 
joining the centers of the ellipse and the circumference is 



r + (a±/2)n=(h/2,L/2). 



(5) 



This leads to the following set of equations 

(L/aj_ - av) 2 [(l - a 2 )v 2 + a 2 } = v 2 (6) 



h/*± = 1 



a 



a L 

v a± 



(7) 



where v = cos (j). Eq.[6]is a quartic equation that must be 
solved together with Eq. [4j It has only one solution for v 
in the range [0, 1]. Once v is known, ft is easily obtained 
from Eq. [7] as a function of the aspect ratio, a. 

The unit cell for this case is shown on the left of 
Fig. ^l). It is a prism with square base of side L/y/2 
and height ft, which contains a couple of particles with 
one of their principal axes (not the axis of revolution) 
along the prism height and the other two parallel to the 
prism base. One particle is placed at one (any) of the 
prism vertices with its axes at the base plane forming 
a 7r/4 angle with any given side of the base. The other 
particle is placed at the prism center, with both principal 



axes parallel to the prism base rotated ir/2 with respect 
to the particle at the vertex. These two particles gen- 
erate the two-folded structure when replicating the cell, 
one yielding layer A and the other layer B, both parallel 
to the prism base. This is shown on the left of Fig. |2]b) 
for oblates with S = y/3. Then, the volume fraction is 
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where we used that the volume of the ellipsoid is v e = 
ircr±cr\\/6. This equation predicts a continuous growth of 
the packing fraction as we increase the asymmetry from 
the value for spheres (S = 1), until it finally reaches a sat- 
uration value for 5 = y/3, where each ellipsoid touches six 
in-plane ellipsoids instead of four. The packing fraction 
obtained at this saturation point is (p m = 0.77073 [1 . 

For S > a/3, Donev et. al. noticed that it is possible 
to stretch the structure obtained for 5 = increasing 
S by the same factor for ellipsoids belonging to different 
layers. Hence, this should lead to any desired aspect 
ratio S > y/3 while preserving a mono-component system 
and keeping the same maximum density. The stretching 
can be done in a direction perpendicular to the layers (in 
general, this would not lead to ellipsoids of revolution), or 
in a direction parallel to any diagonal of the face-centered 
square lattice of layers A or B, since they coincide. In 
the case of our unit cell, this would be in the direction 
of any of the base sides. The relationship between the 
stretching factor, 5, and S is given by [1 
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Note that the stretching leads to an increase of both, 
(7 ram and (J m in- After the stretching, the new (elongated) 
axes, g' and a' in , are now 
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The sides of the prism become L/y/2, sL/y/2, and ft, 
so that its base is now a rectangle of sides L/y/2 and 
sL/y/2. In this case ft can be easily obtained from ft = 
^v e /((p m sL 2 ). The center of both particles are kept in 
the prism vertex and center, and their angles between 
the largest principal axis and the direction parallel to 
the prism side of length L / y/2 turns into 



ip = 1/2 ( 7r =b arctan 



(ii) 



These angles are 7r/4 (minus sign) and 37r/4 (plus sign) 
for s = 1 (no stretching), and tt/2 for s — » oo. The 
stretched cell is shown on the right of Fig. [2] a). In the 
simulations, we set cF min — 1, so the cell and particles are 
rescaled by l/a' min . A simulation snapshot of 120 oblates 
with 5 = 3 is shown on the right of Fig. [5]b). The cells 
before and after the stretching are particular cases of the 
more general sm2 family of structures [U [5] . 
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C. Replica Exchange Monte Carlo 

This technique was developed to enhance sampling at 
difficult (high density / low temperature) conditions [14]- 
[16] . It is based on the definition of an extended en- 
semble whose partition function is given by Qextended — 
riSi Qit being n r the number of ensembles and Qi the 
partition function of ensemble i. This extended ensemble 
is sampled by n r replicas, each replica placed at each en- 
semble. The extended ensemble justifies the introduction 
of swap trial moves between any two replicas, whenever 
the detailed balance condition is satisfied. For hard par- 
ticles, it is convenient to make use of isobaric-isothermal 
ensembles and perform the ensemble expansion in pres- 
sure [17 . For this particular choice, the partition func- 
tion of the extended ensemble turns [T71 H8] 

n r 

Qextended J^J QntPi, (12) 
i=l 

where QattP; is the partition function of the isobaric- 
isothermal ensemble of the system at pressure P^, tem- 
perature T, and with N particles. 

We implemented a standard sampling of the NT Pi en- 
sembles, involving independent trial displacements, ro- 
tations of single ellipsoids, and volume changes. To 
increase the degrees of freedom of our small systems 
(N = 120), we implemented non-orthogonal paral- 
lelepiped cells. Thus, sampling also includes trial changes 
of the angles and relative length sides of the lattice vec- 
tors defining the simulation cell. This is done while 
rescaling the cell sides and particles positions to preserve 
volume and keep a simple acceptation rule. Swap moves 
are performed by setting equal probabilities for choosing 
any adjacent pairs of replicas, and using the following 
acceptance rule [17] 

P acc = min(l,exp[/3(P < - P 3 )(V t - Vj)]), (13) 

where j3 = l/(/c#T) is the reciprocal temperature and 
Vi — Vj is the volume difference between replicas i and j. 
Adjacent pressures should be close enough to provide rea- 
sonable swap acceptance rates between neighboring en- 
sembles. In order to take good advantage of the method, 
the ensemble at the smaller pressure must also ensure 
large jumps in configuration space, so that the higher 
pressure ensembles can be efficiently sampled. 

Simulations started from the Donev-sm2 structures de- 
scribed in the previous section. We first perform about 
5 x 10 12 trial moves at the desired state points, during 
which we check that the replicas have reached a station- 
ary state (thermalization stage). It should be noted that 
achieving a stationary state requires considerably less 
simulation steps when starting from these ordered config- 
urations than when starting from loose random configu- 
ration cells [6]. We then perform 1 x 10 13 additional sam- 
pling trials. Maximum particle displacements, maximum 
rotational displacements, maximum volume changes, and 




FIG. 3. a) Equation of state for oblates with 6 = 1.2, 
Z(ip). b) Isothermal compressibility, xM- c ) Order pa- 
rameter, Qo(cp). The insets zoom in the highlighted regions. 
Arrows point out the fluid-plastic, plastic-sfcc, and sfcc-sm2 
phase transitions. 

maximum changes of the lattice vectors are adjusted for 
each pressure to yield acceptance rates close to 0.3. Since 
an optimal allocation of the replicas should lead to a 
constant swap acceptance rate for all pairs of adjacent 
ensembles [19], we implemented a simple algorithm to 
smoothly adjust the intermediate pressures while keeping 
the maximum and minimum pressures fixed. To start the 
simulations, we use a geometric progression of the pres- 
sure with the replica index. These adjusting procedures 
are performed only during the equilibrating stage. Verlet 
neighbor lists [20] [21] are used to improve performance. 
We set N = 120 ellipsoids and n r — 64 (to cover a wide 
range of densities while keeping large swap acceptance 
rates) . More details are given in previous works [6j [TT] . 

III. RESULTS 

A. Hard ellipsoids phase diagram 

In the initial stage, all simulation replicas start with 
a sm2 structure. Then, we let the system of replicas 
decompress to yield a stationary state. Once this fi- 
nal state is reached, the compressibility factor (Z(cp) = 
pP/p), the isothermal compressibility (x(^) = N(< p 2 > 
— < p > 2 )/ < p > 2 ) and the Qe-order parameter 

(Qe(y) = (if E™=-6 I <Y em (e,4>)> I 2 ) ) are calcu- 
lated. In these expressions p is the particle number den- 
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FIG. 4. Radial distribution functions, g(r), (dark) and their 
corresponding radial orientational order parameter functions, 
p(r), (light gray) for 8 = 1.2, and for Z = 259, 111, and 60, as 
shown in panels a) , b) and c) , respectively. The corresponding 
snapshots are shown in panels a) (sm2), b) (sfcc), and c) 
(plastic crystal) of Fig. [5] respectively. 



sity and <YQ m (0, (j>) > is the average over all bonds and 
configurations of the spherical harmonics of the orien- 
tation polar angles and (j) [7J |22j [23] . Fig. [3] shows all 
these quantities for quasi-spherical oblates (S = 1.2). The 
graphs also include arrows pointing to the phase transi- 
tions. Moving from low to high pressures, we first find a 
fluid-solid phase transition, where an isotropic fluid and 
a plastic crystal coexist [3]. Then, there is a solid-solid 
transition between a plastic crystal and a sfcc crystal. 
Finally, at high density we observe another solid-solid 
transition, between the sfcc structure and the densest 
sm2 crystal. The first two transitions are first order, 
as pointed out by the Z(ip) discontinuity (see Fig. [3] a)) 
and by the formation of bimodal density distributions 
(not shown). The sfcc-sm2 transition is not first order, 
according to the continuous Z(ip) (Fig. [3] a)), the very 
small kink of x(<£>) (Fig. [3] b)) , and the slight deforma- 
tion of the Gaussian density distributions (not shown). 
At this point we must stress that this last conclusion may 
be affected by the small system sizes we are considering. 
The order parameter Qq{^p) also provides evidences of the 
transitions, as shown by the arrows in panel c). Here, a 
very subtle increase of Qe(cp) is observed for the sfcc-sm2 
transition. 

We can get a more clear evidence of the solid-solid 
phase transitions by studying the behavior of the ra- 
dial distributions functions, g(r), and the radial ori- 
entational order parameter functions, p(r), defined as 
p(r) =< (3(uf • u^) 2 - 1) > /2 [24 . They are shown 
in Fig. [4] Panels a), b), and c) of Fig. [4] are built for a 
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FIG. 5. Equilibrium structures for oblates with 5 = 1.2. 
Lower panels show snapshots where ellipsoids belonging to 
the same layer are equally colored. Upper panels show the 
corresponding front views where oblates are represented by 
small plate-like particles to highlight their crystal-like posi- 
tions and orientations. Columns a), b), and c) correspond to 
sm2, sfcc, and plastic solids. Pressure decreases from a) to c). 



high (Z = 259), intermediate (Z = 111), and low pres- 
sure (Z = 60) solids, respectively. Thence, g(r) and p(r) 
of panel a) are ensemble averages of, mostly, sm2 struc- 
tures. For panel b) the dominant structure is a sfcc and 
for panel c) the most frequent is a plastic solid. The 
g(r) for the sm2 structure shows two main peaks, a first 
one at a distance close but larger than cF mini and a sec- 
ond larger one, at a somewhat larger distance. These 
two peaks correspond to the first coordination shell, for 
the intra (four) and extra-plane (eight) neighbors, respec- 
tively (see Fig. [I]). As shown by the p(r) function, the 
first p(r) ~ 1 and the second p(r) ~ —0.5 peaks cor- 
respond to the parallel (intra-plane) and perpendicular 
(extra-plane) orientations, respectively. Thus, in general 
for this arrangement, a positive p(r) is associated with 
g(r) peaks for par ticles belonging to the same plane (A 
or B) (see section II B| ), whereas a negative p(r) corre- 
sponds to correlations between particles of planes A and 
B. Contrasting with the p(r) function of Fig. [4] a), panel 

b) shows a positive p(r) for all distances, r. This implies 
a background parallel long-range orientational particle- 
particle correlation. Furthermore, the distance between 
the first and second g(r) peaks enlarges, pointing out the 
differentiation between the stretched and the unstretched 
sides of the sfcc cell. Finally, panel c) shows the g(r) and 
p(r) functions of a plastic crystal. That is, while the g(r) 
still shows the well-defined peaks of a solid (fcc-like), p(r) 
shows no long-range angular correlations. 

All features of the g(r) and p(r) functions of Fig. [4] 
correspond to the structures shown in Fig. [5] The lower 
panels of Fig. [5] illustrate snapshots of the sm2, sfcc, and 
plastic solids, in correspondence with panels a), b), and 

c) of Fig. [4] To gain clarity, the upper panels of this figure 
show the corresponding front views where the oblates are 
represented by small plate-like particles. Also, particles 
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FIG. 6. Phase diagram of hard ellipsoids of revolution. 
The spherical case is given for (5=1, whereas prolates are 
at the left and oblates at the right. The dark solid line 
is the maximally achievable density pQ. There are several 
transition types. These are: Isotropic-nematic fluid-fluid (as- 
terisks), isotropic-plastic fluid-solid (squares), nematic-sm2 
fluid-solid (diamonds), isotropic-sm2 fluid-solid (upward tri- 
angles), plastic-sm2 solid-solid (downward triangles), plastic- 
sfcc solid-solid (circles), and sfcc-sm2 solid-solid (crosses). 
Pairs of open and solid symbols are employed to show co- 
existence regions. Single symbols are employed to point out 
higher order transitions. The inset zooms in the sfcc stable 
region (in-between solid circles and crosses). 



in the lower panels are colored according to their posi- 
tions, to highlight the layering of the structures. From 
this pictures, the sm2, sffc, and plastic-solid structures 
can be easily recognized. 

The existence of a subtle sfcc-sm2 transition for oblates 
with 5 = 1.2 encouraged us to explore the boundaries of 
the sfcc stable region. For this purpose, we performed 
a similar analysis, but now considering other aspect ra- 
dios for both oblate and prolate cases. The obtained re- 
sults allow the construction of a refreshed hard-ellipsoid 
phase diagram, which is shown in Fig. [6j The diagram 
is split in half by the hard sphere case (6 = 1). Prolate 
cases are at the left and oblate systems are at the right 
of this vertical line. Particles' asymmetry increases by 
moving away from this central line. The extreme cases 
are infinitely narrow needles (1/5 — » at the left) and 
infinitely thin plates (1/5 —> at the right). Prolates 
with 5 > 3 are not included since a larger number of par- 
ticles is needed to fulfill the minimum-image convention 
for all densities. At high densities, we are placing the 
(currently accepted) maximally achievable density pQ as 
a black solid line. All transitions found in our simula- 
tions are indicated in this chart. Isotropic-nematic fluid- 
fluid transitions are given as asterisks, isotropic-plastic 
fluid-solid transitions as squares, nematic-sm2 fluid-solid 
transitions as diamonds, isotropic-sm2 fluid-solid transi- 
tions as upward triangles, plastic-sm2 solid-solid transi- 
tions as downward triangles, plastic-sfcc solid-solid tran- 



FIG. 7. Pressures at which transitions take place as 
a function of £ _1 for prolates (left) and oblates (right). 
Different symbols match the transitions given in Fig. [6] 
These are: Isotropic-nematic fluid- fluid (asterisks), isotropic- 
plastic fluid-solid (squares), nematic-sm2 fluid-solid (dia- 
monds), isotropic-sm2 fluid-solid (upward triangles), plastic- 
sm2 solid-solid (downward triangles), plastic-sfcc solid-solid 
(circles), and sfcc-sm2 solid-solid (crosses). 



sitions as circles, and sfcc-sm2 solid-solid transitions as 
crosses. Isotropic-nematic fluid-fluid and sfcc-sm2 solid- 
solid transitions are higher order, and thus, they are 
shown as a single point, located at the packing fraction 
where the isothermal compressibility reaches a local peak. 
All other transitions are first order, and thus, a couple 
of symbols are used to denote the borders of the coexis- 
tence regions. The phase boundaries are determined by 
means of the histogram re- weighting technique described 
elsewhere [25 , 26 . The numerical values are given in 
tables [T] and [TTJ According to previous results, increas- 
ing the system size would slightly shift the transitions 
towards larger densities and would make them slightly 
wider [IU Hi]. 

The inset of Fig. [6] illustrates a zoom of the high den- 
sity area above the plastic region. There, it can be ap- 
preciated a narrow density region where the sfcc solid 
spontaneously forms. Hence, for very low asymmetries, 
the plastic solid turns into a sfcc solid before taking the 
form of a sm2 structure under very high compression. 
It should be noted that the sfcc solid region is not so 
small in terms of the pressure range at which it is sta- 
ble. As shown in Fig. [7| the sfcc stable region can be 
several hundreds of /3Po~^ n wide (see also table [i]). The 
existence of a stable sfcc structure above the plastic solid 
region is in agreement with Radu et. al. [5 . Nonethe- 
less, and despite that the phase boundaries for the sfcc 
structure were not given, it was suggested that this struc- 
ture should be the most stable for 1/5 < 0.65. Our data 
show that this is only true for a very small density re- 
gion above the plastic solid, so a sfcc-sm2 transition ap- 
pears at very high pressures and densities. Free energy 
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TABLE I. Coexistence volume fraction borders ip and pressure P* = /3Pcr^ n of transitions for cases with low asymmetry. 
Subindexes / and h refer to the low and high density borders. A dash means the absence of the transition. Errors are always 
below 3% for all quantities. 
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plastic-sm2 


plastic-sfcc 


sfcc-sm2 


a 
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p* 


(pi 
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p* 


(f 


p* 


1.400 


0.572 


0.587 


18.3 








0.656 


0.666 


47.4 


0.675 


48.9 


1.300 


0.532 


0.555 


12.9 








0.679 


0.687 


69.7 


0.699 


85.6 


1.200 


0.508 


0.548 


11.5 








0.699 


0.704 


116 


0.717 


172 


1.150 


0.501 


0.542 


10.9 


_ 


_ 


_ 


0.708 


0.713 


156 


0.725 


269 


1.100 


0.495 


0.539 


10.8 








0.720 


0.723 


265 


0.732 


535 


1.050 


0.492 


0.538 


10.8 








0.731 


0.733 


634 


0.737 


1528 


1.000 


0.490 


0.537 


10.5 








0.740 


0.740 


oo 


0.740 


CO 


0.952 


0.491 


0.537 


10.3 








0.732 


0.733 


661 


0.739 


2412 


0.909 


0.495 


0.540 


9.91 








0.721 


0.724 


256 


0.734 


605 


0.870 


0.499 


0.540 


9.40 








0.711 


0.716 


152 


0.726 


251 


0.833 


0.506 


0.541 


9.08 








0.701 


0.706 


100 


0.717 


141 


0.769 


0.527 


0.552 


9.48 








0.679 


0.688 


54.3 


0.694 


59.0 


0.714 


0.561 


0.572 


11.2 


0.657 


0.674 


33.1 












0.667 


0.614 


0.614 


16.3 


0.638 


0.651 


20.2 













TABLE II. Coexistence volume fraction borders if and pressure P* = /3Pcr^ n of transitions for cases with large asymmetry. 
Subindexes / and h refer to the low and high density borders. A dash means the absence of the transition whereas a blank 
means that the experiment was not carried out. Errors are always below 4% for all quantities. 
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0.384 
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4.000 


0.452 


3.14 














3.000 


0.543 


8.30 








0.566 


0.618 


9.50 


2.000 






0.579 


0.625 


15.3 








1.732 






0.595 


0.637 


19.3 








1.500 






0.648 


0.665 


35.4 








0.577 






0.589 


0.632 


10.5 








0.500 


0.561 


6.52 








0.565 


0.612 


6.91 


0.333 


0.499 


1.87 








0.566 


0.609 


2.99 


0.250 


0.412 


0.550 








0.565 


0.606 


1.61 


0.200 


0.348 


0.229 








0.577 


0.612 


1.09 



calculations through thermodynamic integration should 
confirm this finding [27-29 . As mentioned in the intro- 
duction, this transition should necessarily exist in order 
to be consistent with the fact that the sm2 structure is 
the one leading to the maximally achievable density for 
all aspect ratios. 



B. Comparing runs from random and ordered 
initial conditions 

This section is devoted to compare the REMC results 
obtained by starting from loose random cells (data taken 
from [6] ) with those obtained from dense sm2 initial con- 
figurations. In order to support the ergodicity of the 
simulated systems, one should obtain the same results re- 
gardless of the initial configurations. For certain systems 
at very high pressures, the ergodic condition is difficult to 
achieve. This is simply because simulations tend to get 
stuck on configuration space at high densities. When this 



occurs, in turn, dynamical properties usually show solid- 
like behaviors. That is, particles relaxation times diverge, 
and diffusion coefficients turn practically zero [30]. Thus, 
with the aim of enhancing sampling at these difficult con- 
ditions, certain Monte Carlo techniques were developed, 
usually introducing unnatural displacements while keep- 
ing the detailed balance condition. This is the case of the 
REMC technique, where replicas can travel throughout 
different ensembles. 

As mentioned in section |II C[ an improved sampling 
can be obtained whenever a relatively large swap accep- 
tance rate is achieved between all set pressures. Nonethe- 
less, this does not guaranty ergodicity. For instance, 
replicas may frequently jump from one side to the other 
of a coexistence region without residing during a large 
number of steps where they "do not belong". That is, 
the replica having the most compressed fluid-like config- 
uration may swap positions with the replica having the 
less compressed solid-like structure, across a coexistence 
region, but this fluctuation may not last long enough to 
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FIG. 8. Oblates with 5 = 1.3. a) Compressibility factor 
Z(ip). b) Isothermal compressibility, xM- c ) Bond order pa- 
rameter, Qq(lp). Dark circles correspond to simulations start- 
ing from dense sm2 structures, whereas light squares corre- 
spond to simulations starting from loose random cells (data 
taken from [6]). 



FIG. 9. Prolates with 8 — 1.5. a) Compressibility fac- 
tor Z(ip). b) Isothermal compressibility, x(^). c) Bond or- 
der parameter, Qo(cp). Dark circles correspond to simulations 
starting from dense sm2 structures, whereas light squares cor- 
respond to simulations starting from loose random cells (data 
taken from [6]). 



allow the fluid-like configuration become solid-like and 
vice-versa. Consequently, there would be a low rate of 
generation of new solid-like structures, and the improved 
sampling may be not sufficient to attain ergodic condi- 
tions. 

An example where REMC works well is given in Fig. [8] 
There it is shown how the results obtained starting from 
dense sm2 structures (dark circle curves) match the ones 
started from loose random initial conditions (light square 
curves). This figure corresponds to oblates with S — 1.3, 
and considers the pressure range 2.0 < /?Pcr^ n < 200 
for random initial conditions and 8.5 < /3Pa^ in < 1000 
for sm2 initial conditions. As observed, both runs present 
the same density jumps, pointing out the isotropic-plastic 
fluid-solid, and the plastic-sfcc solid-solid transitions. In 
fact, the dark x(y?) curve also shows a sfcc-sm2 transi- 
tion, appearing as a small shoulder developed at the right 
of the plastic-sfcc transition peak. This subtle transi- 
tion, which is better seen from the slight distortion of the 
Gaussian density distributions and from the radial distri- 
bution functions (not shown), is not captured by the run 
with random initial conditions. Nevertheless, the fact 
that the plastic-sfcc solid-solid transition is captured by 
both runs is quite remarkable, given the extremely high 
density at which it occurs. In general, we observed good 
agreements between runs started from different condi- 
tions for ellipsoids with 5 < 1.5. 

The comparison for prolates with 5 = 1.5 is shown 
in Fig. [9] Here, it is seen a good match between the 
two cases only for densities below the isotropic-sm2 fluid- 



solid transition. In fact, the transition is not captured 
when starting from random initial configurations. In- 
stead, the fluid curve extends towards higher densities, 
defining a metastable branch which cannot be broken 
even by a very large number of REMC cycles. This 
metastable branch is analogous but stronger than the one 
for hard spheres [3TH33] . It should be emphasized that 
the relationship X = N((p 2 ) - (p) 2 )/(p) 2 = dp/d(/3P) 
is satisfied for both data sets. Fulfilling the above rela- 
tionship, suggested to hold only at equilibrium [34], is 
then a necessary but not sufficient condition for equi- 
librium [35]. Note that the order parameter Qq indi- 
cates the development of some bond order at the tran- 
sition. Another curiosity is that an extrapolation of the 
metastable branch would intersect the solid branch. It 
should be said that for densities below but close to the 
transition, the isotropic fluid shows extremely low trans- 
lational and rotational diffusion coefficients [36-38 , i. e. 
dynamics turns glassy. Thus, when starting from random 
initial configurations, the glassy dynamics not only hin- 
ders the solid formation when dynamics simulations are 
used, but also when REMC is the employed technique. 
The other way around, when REMC produces an ap- 
parent stationary state which differs from equilibrium at 
high densities, a glassy dynamics may be expected. The 
occurrence of disordered structures at very large densi- 
ties was also found experimentally [39 and by means of 
computations [7J [40] . These references report that disor- 
dered states yield a maximally random jammed density, 
if « 0.712, peaking at 5 = 1.5. Additionally, they show 
that prolates produce slightly higher maximally random 
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FIG. 10. Oblates with 6 = 4. a) Compressibility factor Z((p). 
b) Isothermal compressibility, xi}P)- c ) Bond order parame- 
ter, Qq(lp). Dark circles correspond to simulations starting 
from dense sm2 structures, whereas light squares correspond 
to simulations starting from lose random cells (data taken 
from j6]). 



jammed densities than oblates. 

Another example where curves do not entirely agree is 
shown in Fig. [lOj Here again, a good match between the 
different runs is obtained for densities below the nematic- 
sm2 fluid-solid transition (the agreement on the descrip- 
tion of the isotropic-nematic fluid-fluid transition is very 
good). For densities above the fluid-solid transition and 
for a given pressure, the run with initial random con- 
figurations yields smaller densities, isothermal compress- 
ibilities, and Qq values than the ones obtained starting 
from sm2 structures. Nonetheless, both runs capture a 
fluid-solid transition and, in addition, the resulting solid 
structures are similar. We then justify differences due to 
the appearance of imperfect sm2 structures when starting 
from disordered configurations. This turns evident from 
a snapshot overview (not shown). It should be mentioned 
that only certain box-shapes can hold an exact number 
of particles of a given perfect sm2 structure. This ef- 
fect, for the system sizes we are employing, is important. 
In previous work [6] we only let free the angles of the 
simulation cells. However, this procedure seems to be 
not enough to completely relax the solid phase. In the 
present work we are letting angles and sides to vary in- 
dependently. Hence, systems with different degrees of 
freedom are being compared and this is why differences 
appear for the solid phases. 

Summing up, starting from random initial configu- 
rations has the following characteristics: a) It natu- 
rally produces highly dense structures which may assist 
proposing the equilibrium structure or confirm/refute the 



believed equilibrium structure, b) It allows the detec- 
tion of strong metastable regions, which may be linked 
to systems showing extremely low dynamics. Thus, when 
these metastable regions appear for dense fluid-like struc- 
tures, a glassy behavior can be expected, c) Long runs 
are frequently necessary to reach a stationary state. On 
the other hand, starting from a-priori set dense crystal 
structure has the following characteristics: a) When the 
crystal structure corresponds to equilibrium at high den- 
sities, the obtained high density branch should also cor- 
respond to equilibrium, while transitions are easily cap- 
tured by decompression, b) Thus, relatively short runs 
would produce equilibrium since metastable regions are 
avoided, c) The imposition of a crystal structure not 
representative of equilibrium may lead to a stationary 
state inside a strong metastable region (strategies to find 
structure candidates to reach maximally achievable den- 
sities, i. e. structure candidates for equilibrium at very 
high pressures, are given elsewhere [41-43 ). From com- 
bining random and crystal initial conditions a better un- 
derstanding of the system is possible, i. e. by comparing 
results from both strategies one can detect metastable 
regions or support ergodicity. 



IV. CONCLUSIONS 

Further details on the phase diagram of hard ellipsoids 
of revolution were captured by decompressing high den- 
sity sm2 structures by means of replica exchange Monte 
Carlo simulations. Mainly, we observed a very rich be- 
havior for quasi- spherical oblates and prolates. These 
systems, from low to high pressures, show the following 
phases: isotropic fluid, plastic solid, stretched-fcc solid, 
and sm2 solid. According to our data obtained for small 
system sizes, the first three transitions are first order, 
whereas the last one is a subtle, high order transition. 
This picture is consistent with the fact of having the sm2 
structure capable of producing the maximally achievable 
density. 

Replica exchange Monte Carlo (REMC) simulations 
started from dense sm2 structures produce, by decom- 
pression of the initial configurations, equilibrium states 
for all set pressures. This, in addition, is yield in rel- 
atively few Monte Carlo steps. Thus, when there are 
structure candidates for the maximally achievable den- 
sity, REMC simulations started by setting them as start- 
ing cells would provide a fast way to produce the whole 
phase diagram. 
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